Retention capacity of random surfaces 
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We introduce a "water retention" model for liquids captured on a random surface with open 
boundaries, and investigate the model for both continuous and discrete surface heights 0, 1, ... n — 1 
on a square lattice with a square boundary. The model is found to have several intriguing features, 
including a nonmonotonic dependence of the retention on the number of levels: for many n, the 
retention is counterintuitively greater than that of an n + 1-level system. The behavior is explained 
using percolation theory, by mapping it to a 2-level system with variable probability. Results in one 
dimension are also found. 
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Consider a bounded horizontal random surface with a 
landscape of varying height, as shown in Fig.[l] A liquid 
such as water is dripped over the surface and is allowed 
to drain out all of the boundaries. Internal sites in val- 
leys capture the water and create ponds, and eventually 
all the ponds fill up to their maximum height. We are 
interested in finding the total amount of water retained 
in the system when the maximum heights are reached. 
Physically, this problem is related to coatings on a ran- 
dom surface and the properties of landscapes and wa- 
tersheds. Theoretically, it is related to the topology of 
random surfaces [H |2] and to invasion percolation (IP), 
but with some interesting new features. 

We study this problem on a regular square lattice with 
random heights assigned to each site. The systems are 
square of size L x L with draining boundaries on all four 
sides. Extensive simulations were performed with uni- 
formly distributed discrete heights 0, l,2,...n — 1 for 
values of n ranging from 2 to 100, and also for a con- 
tinuum of heights 0, 1. We also studied a 2-level system 
with variable occupation probabilities of the 2 heights. 
The simulation method we used is a form of IP in which 
we effectively reversed the flow and flooded the system 
from the outside with higher water levels, and recorded 
the level of the water in a pond when it was first flooded. 
The retention is the difference between that level and the 
height of the terrain below the pond. 

Fig. [2] shows the average retention i?n^^ on n-level sys- 
tems, for n = 2, ... 8, as a function of L. Here we assume 
that all terrain heights occur with equal probability. As 
expected, the retention grows as for large L, and gen- 
erally grows with n, as more levels create deeper ponds. 
However, we found deviations to this expected behavior. 
As seen in Fig. [2j there is a crossover in the curves for 



n = 2 and 3: for small L, i?^^^ < but for L > 51, 



i?2^'^ > i^3^^ This is in spite of the fact that 3- level sys- 
tems can have ponds of water of depth 2, while 2-level 



systems can only have ponds of depth 1. Further study 
shows additional crossings between levels n and n + 1 at 
seemingly random n's, and at larger values of L (Table 

0. 



II DlBc oD3Bgl 





FIG. 1. (color online) Examples of water retention on square 
systems with an equal distribution of terrain heights. Top left: 
heights through 9; Top right: heights and 1. Green (gray 
with black height numbers) are dry and above water, blue 
(gray with withe numbers) are filled with water, and yellow 
(light shade) are spillover sites of the height of an adjacent 
pond (and only shown if neighboring a wet site). Bottom: 
Perspective view of a water- filled 10-level system, with dry 
sites shown as grey (lighter shading on top). 

In this Letter we explain some of the puzzling features 
of this model, though many questions remain. Some re- 
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lated issues, especially involving multi-level nonuniform 
systems, are discussed in [3 . 

To analyze the multi-level discrete model, we make a 
decomposition of in terms of the retention in a 2- 
level system with varying p, R^\p)^ where p = Prob(O) 
is the probability or fraction of sites with terrain height 
in the 2-level system: 

n-l 
i=l 

The i = 1 term represents the amount of water retained 
up to just the first level, for which all sites of terrain 
height 1 or higher can be considered as level 1. The 
net fraction of 0- height sites is 1/n. The i = 2 term 
represents the total amount of water at just the second 
level; here we collapse the first 2 levels into the new level 
(fraction 2/n), and the rest of the sites can be considered 
as level 1. Likewise, the remaining terms follow. 
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FIG. 2. (color online) Retention l6a^ vs. L for n = 2, 3, . . . , 8 
with a uniform distribution of levels, showing the crossing of 
the curves for n — 2 and 3 at L 51. Additional crossings 
are found for larger L (see Table [l]) . For small L there is an 
ordering of the retention with n, which is violated for larger 
L. 

It is also possible to show that {l/i)R^\i/n) is equal 
to the number of sites with retention i. Thus, the total 
number of wet sites is 

n— 1 -J 
i=l ^ 

It is therefore necessary to just know the behavior of 
the 2-level system with a varying p in order to predict 
the behavior of all n-level systems (including those with 
nonuniform level distributions). We have carried out sim- 
ulation of R^2^ {p) for p = 0.01, 0.02, . . . , 0.99 for various 
L, and the results are shown in Fig. |3] For small L, the 
curve is rounded and peaked close to P(0) = 1/2, but 
for larger L it approaches a ramp of slope 1, up to the 



value of p = Pc = 0.592 746 (the site percolation thresh- 
old), after which it drops off precipitously and rapidly 
approaches 0. This behavior is best understood from the 

limit L ^ oo. Define r2{p) = lim/^^oo R^\p)/ L'^ • Iii an 
infinite system, all finite clusters of 0-sites retain water 
and the infinite cluster alone drains off. Thus, the total 
retention is 

r2(p)=p-Poo(p), (3) 

where Poo{p) is the fraction of sites belonging to the 
infinite cluster. Very close to (and above) 

Pc 1 Poo ^ 

a{p — PcY + • • • where a is a constant and {3 = 5/36 
[4 . Because Poo rises quickly as p increases beyond Pc, 
we get the quick drop to zero in r2{p). For finite L, fi- 
nite clusters at the boundaries drain as well, yielding the 
observed finite-size effects. 

Exactly at Pc^ the drainage area is fractal, yet the re- 
tained water is still proportional to for large L with 
corrections proportional to L^^ where df = 91/48 is the 
fractal dimension. We verified that at Pc^ the size dis- 
tribution of the draining clusters (boundary clusters in 
percolation) satisfies n'^ ^ with = 1/df — 1 = 
139/91 ^ 1.527 as predicted in [5 . Our measurement of 
r' = 1.5256 ± 0.003 confirms this prediction to about 20 
times the precision given in [5]. 
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FIG. 3. (color online) Retention fraction R^^\p)/L^ in a 2- 
level system with varying p =Prob(0), for L ranging from 50 
to 1000, averaged over 10 000 to 1000 000 samples for each 
point. 

As a good approximation for large L, we can ignore the 
small contribution to R^^^ (p) for p > Pc and approximate 
R^^\p) = pL^ for p < Pc. Then, from (1), we find the 
following formula for the n-level retention in the large-L 
limit: 

(4) 

where n* is the largest integer such that n*/n is less than 
Pc. Thus, we have i?^^^ - (l/2)P^ and R^^^ - (l/3)P^ 
which indeed gives R^^^ > Rf^ for large L. This result 
can be explained simply by the fact that for the 2-level 
system, roughly half the sites are O's and filled with wa- 
ter, while for the 3- level system, only 1/3 of the sites 
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are O's; very few ponds are filled to a level of 2 because 
those sites correspond to clusters above the percolation 
threshold. 

To explain the crossing, we must also explain why the 
curves for R^"^ are ordered i^2^^ < R^^^ < R^\.. for 
small L. This can be understood qualitatively from the 
behavior of R^\p) for small L as in Fig. bl because 
those curves are smooth, equation ([T]) will be a gradual, 
increasing function of n. To be more rigorous, we con- 
sider the smallest system possible: a 3 x 3 system, which 
has only one site that can hold water (the center site), 
and only four sites that can block it, as the corner sites 
are irrelevant. A direct calculation yields 

R^n^ = {ri^ - l)(3n' - 2)/(60n^) , (5) 

which is a monotonically increasing function of n. (De- 
tails of the derivation will be given in a future paper.) 
Because the ordering is verified for L = 3 but not for 
large L, crossing must necessarily occur for some L. 

The curves that are "out of order" and cross are those 
in which = n^{n^ + l)/(2n) is greater than r^+i, 
by (|4|. This occurs when the fractional part of Pcfi is 
between and 1 — Pc ^ 0.407. The crossing curves 
(n,n + l) are at (2,3), (4,5), (7,8), (9,10), (12,13), 
(14, 15), (17, 18), etc. We have verified the first six cross- 
ings as shown in Table |lj For n > 30, the simple analysis 
based upon Q evidently breaks down as contributions 
from R2{p) for p > Pc become important, and the cross- 
ings are predicted to become less frequent, though we 
have not measured them directly. 

TABLE I. Crossing points where Ri^ — R^^+i^ extrapolated 
to non-integer L*. 



n and (n + 1) 


L* 




2 and 3 


51.12 


790 


4 and 5 


198.1 


26 000 


7 and 8 


440.3 


246 300 


9 and 10 


559.1 


502 000 


12 and 13 


1390.6 


4 288 500 


14 and 15 


1016.3 


2 607 000 



There are also small ponds with higher levels, corre- 
sponding to clusters above the threshold; these allow the 
average water level to be above Pc- Fig. [4] shows the wa- 
ter level of sites when first fiooded as a function of the 
number of sites fiooded, showing the small contribution 
of the ponds of high level. 

In fact, taking the continuum limit of ([T]) we can 
calculate the total retention r per site in the contin- 
uum system directly by integrating the curve of r2(p), 
r = Jq r2{p)dp. The triangular part below Pc gives 
pI/2 exactly, and the tail above Pc gives a small cor- 
rection. The tail's area extrapolates to 0.0063 for large 
L, and predicts r = p^/2 + 0.0063 = 0.1820, which 
we verified directly to ±0.0002 by measuring the reten- 
tion for systems of up to L = 12 000 and extrapolat- 
ing to L = oo. The retention per site r is equal to 
{wh)/2 ^ wh/2 ^ where w = W^^^/L^. Note, 

^ ~ lo i'^'^ip)/p)^P follows from ^ in the continuum 
limit, and we find w = Pc-\-0.0100 = 0.6028 in agreement 
with direct measurement (see h = 2r/w — 0.6039 

was also independently measured. 

In fact, applying ([3| to Eqs. ([2| and ([T]) we see that 
w — 1 — r = 1/2 — /io 7 and more generally, we find 
the moments of the retention as 

where /ig = x^Poo{x) dx is the q-th moment of Poo- 
Thus, we have found that the moments of Poo assume 
a specific physical interpretation in the context of the 
retention problem. 
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In the limit that the number of levels becomes infinite, 
the discrete system goes over to the continuum one. Now, 
as in traditional IP 0, the fiuid fiows over the lowest 
barrier site on the perimeter of a pond. For a continuum 
bond IP system, the "raining" IP problem has recently 
been considered in [TllH], and the pond-size distribution, 
away from the boundaries, was investigated. 

Here, considering the continuum site system, we find 
that water rises to an average height of h ^ 0.6039 (av- 
eraged over wet sites only, for P ^ oo), which is slightly 
above Pc. The large ponds have a water level that is 
slightly below Pc^ because higher levels produce large 
percolation clusters that would run into the boundary. 



FIG. 4. Level when sites are first flooded vs. the fraction of 
sites flooded on a 1000 x 1000 continuum system. 

The asymptotic behavior of W'^^^ and P^^^ in the con- 
tinuum system is found to be 

W^^^ - 0.6028P^ - 2.4P^-2^ 
P^^) - 0.1820P^ - 1.27P^-^^ (7) 

where the second terms reflect the effects of the drainage 
sites and ponds of lower water level near the boundary. 
At each level in a discrete system, the drainage area is 
just all the clusters touching the boundary, which ex- 
tends into the system a distance of the correlation length 
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^ ~ \p — Pc\~^' Integrating this over p up to Pc — cL~^^^ 
and multiplying by the perimeter 4L gives a depletion 
zone cx = L^/^. The value of the exponent 1.25 

was verified numerically to ±0.05. Recently, it has been 
shown that watersheds, bounded by the "continental di- 
vide" between drainage regions, have a similar fractal 
dimension df ^ 1.22 j^. It appears that these two prob- 
lems, however, are different, despite the similarity of ex- 
ponents. 

The assignment of a terrain height for each site us- 
ing a probability corresponds to a grand canonical type 
of description. One can also distribute the levels canoni- 
cally, with exactly 1/n of them of each height. We carried 
out simulations using this ensemble and found only small 
differences. For L = 3, 4, and 5, we also carried out an 
exact enumeration of all canonical states. For small L, 
A = [canonical] — [grandcanonical] > 0, and for 
larger L, A decreases. For n = 2, A appears to approach 
0, while for n = 3, 4, 5, A appears to approach a negative 
constant for large L. Because the value of the retention 
itself grows as L^, the relative difference A/i?l^^ is very 
small. We verified that changing to the canonical en- 
semble does not affect the crossing behavior of the R^^^ 
curves. 

We studied the distribution rig of ponds of s sites and 
verified that the system self-organizes to the percolation 
critical point with Ug ~ and r = 187/91. Unlike 
standard percolation, we cannot write exact formulas for 
any — not even for ni. However, we can make an esti- 
mate for ni as follows: The probability that a site is in 
a pond of size 1, of water height between x and x -\- dx^ 
is given approximately by 

Pi{x)dx = 4(1 - xfx[l - (1 - xf]dx (8) 

where the factor of 4(1 — x)^ is the probability that 3 of 
the neighbors are of higher terrain height (4 possibilities), 
X is the probability that the site itself is of terrain height 
less than or equal to x, and [1 — {1 — x)^] is the probability 
that the spillway site has at least one neighbor lower than 
X, so the spillway can drain at least to the next sites. This 
gives ni ^ /J^ Pi{x)dx = 0.01624 where xq = h = 0.6039 
is the average water height surrounding the cluster. This 
compares to a measurement of rii = 0.015 95. Like- 
wise, the average water height of the ponds of size 1, 
hi = /J^ xPi{x)dx/ J^^ Pi{x)dx = 0.6904, is close to the 
extrapolated measured value 0.6887. 

We also studied the model in one dimension (ID), 
where there are no crossings, however exact results for 
all quantities can be found. Consider, for example, the 
semi-infinite line (sites 1,2,3,...) so that water can spill 
only through the left edge, and assume a uniform distri- 
bution of barriers in [0, 1]. As we look at the ID ponds, 
starting from the left edge, the water level keeps rising 
the farther we venture into the line. In fact, each pond 
begins when a record-height barrier is encountered, and 



ends when the next, yet higher barrier, is met. 

The probability that the barrier at site /c + 1 is taller 
than all the k preceding barriers is x^ dx = 1/ {k -\- 1). 
This is also the probability that a pond starts (or ends) 
at site k-\-l. Because the barriers demarcating the ponds 
occur with probability l//c at site /c, it follows that the 
typical size of ponds, k sites away from the edge, is k. The 
ponds grow linearly with their distance from the edge. 

Next consider Ps{k)^ the probability that a ID pond of 
size s is k sites away from the edge, in sites k -\- l^k -\- 
2, . . . , /c + s. For that pond to have water level x, the first 
k — 1 sites must have barriers lower than x (with proba- 
bility x^~^), as do sites /c + l,/c + 2, + 5 (probability 
x^). Site k contains a barrier of height x (probability 
dx), and site k-\- s-\-l contains a barrier of height y > x 
(probability 1 — x). Thus the probability for a pond of 
level X is x^~^+^(l — x) dx. Integrating over x, we obtain 
the required probability: 

Note that ^^iPs{k) = l/(/c + 1), consistent with our 
previous result, and that the moments of Ps{k) diverge, 
which is why we estimated the typical pond size instead. 

Similarly, the probability for having k draining sites at 
the edge is 

P,rUk) = / = (fcTI)! ■ 

These results illuminate the analogous quantities in 2D, 
where however no exact results could be found. 

In this Letter we have only touched upon the questions 
that one may ask about the retention model. There are 
many more questions that are unsolved, including the 
exact results for the size distribution of the clusters, the 
average retention as a function of the distance from an 
edge, the behavior on other lattices, on systems with dif- 
ferent boundary shapes, in higher dimensions, and sys- 
tems with a tilt. We believe it is an interesting model 
that warrants much further study. 

We mention finally that the water retention problem 
was previously studied in the context of surfaces created 
by magic squares [10]. The application to random sur- 
faces is an example of the deeper connections of this prob- 
lem. 
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Note added: This paper was published in Physical Review Letters 108, 045703 (January 25, 2012). Some additional 
material is given in the web page 
Ihttp : //en . wikipedia . org/wiki/Water_retention_on_inatheinatical_surf aces' 

where retention on magic squares is also discussed. Measurements of the crossing points Ri^ ^ = rI^^I for several 
larger values of n have been found and will be given in a future publication. 
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